clear
set mem 80m
set more off
set matsize 800
capture log close

cd "D:\MyPapers\Happiness\TWINSDATA\code\2023May-Dataverse"
use twins_PNAS.dta, clear

log using TwinsIncomeHappinessMainAnalysis.txt, replace

global citylist "d_Chengdu d_Chongqing d_Haerbin d_Hefei"

global alllist "citycode d_Chengdu d_Chongqing d_Haerbin d_Hefei happy_4s twcode mz dz happiness happy_d income income_num income_bracket IVincome age agesquared male BirthOrder_ElderDummy birth_weight VeryEarly_disease_d early_disease_d edu working work_over_40h  married divorced widowed population birth_year health edu_spouse edu_spouse_atmarriage income_spouse_atmarriage Dhappy_4s Dhappiness Dhappy_d Dincome IVDincome sDincome sIVDincome DBirthOrder_ElderDummy Dbirth_weight DVeryEarly_disease_d Dearly_disease_d Dedu Dworking Dwork_over_40h Dmarried Ddivorced Dwidowed Dpopulation Dhealth Dedu_spouse Dedu_spouse_atmarriage Dincome_spouse_atmarriage"

global reshape_list "citycode d_Chengdu d_Chongqing d_Haerbin d_Hefei happy_4s happiness happy_d income income_num income_bracket IVincome age agesquared male BirthOrder_ElderDummy birth_weight VeryEarly_disease_d early_disease_d edu working work_over_40h  married divorced widowed population birth_year health edu_spouse edu_spouse_atmarriage income_spouse_atmarriage Dhappy_4s Dhappiness Dhappy_d Dincome IVDincome sDincome sIVDincome DBirthOrder_ElderDummy Dbirth_weight DVeryEarly_disease_d Dearly_disease_d Dedu Dworking Dwork_over_40h Dmarried Ddivorced Dwidowed Dpopulation Dhealth Dedu_spouse Dedu_spouse_atmarriage Dincome_spouse_atmarriage"

//Correlation comparisons in Section "Choice of IV" 
quiet reshape wide $reshape_list, i(twinpair) j(twcode)
pwcorr income1 IVincome1 income2 IVincome2
sum   income1 IVincome1 income2 IVincome2
quiet reshape long

//Tables S6 and S7: whether within-twin variance of income and happiness are sufficiently large
quiet reshape wide

tab happiness1 happiness2 if mz
tab happy_4s1 happy_4s2 if mz
tab income_bracket1 income_bracket2 if mz

tab happiness1 happiness2 if dz
tab happy_4s1 happy_4s2 if dz
tab income_bracket1 income_bracket2 if dz

quiet reshape long

bys mz: count if abs(Dincome)>=0.2
bys mz: count

//SI S1 "Alternative IVFE Specification": Correlation of First Stage of IVFE-2
pwcorr sDincome sIVDincome if mz, sig star(.05)  

//SI S1 "Alternative IVFE Specification": evidence of systematic under-reporting of co-twin sibling income
ttest income=IVincome

ttest income=IVincome if twcode==1
ttest income=IVincome if twcode==2

ttest income=IVincome if mz
ttest income=IVincome if dz

ttest Dincome=0
ttest IVDincome=0
ttest sDincome=0
ttest sIVDincome=0

//Figures
//SI, Figure S1: Happy vs. Numberic Income
egen bins=cut(income_num),at(0(200)30000)
replace bins=bins+100

bysort bins: egen happy_4s_mean=mean(happy_4s)

label variable happy_4s_mean "Mean Happiness"
twoway (scatter happy_4s_mean bins) (lfit happy_4s income_num)(qfit happy_4s income_num), ///
       ytitle(Mean Happiness) xtitle(Income) saving(MeanHappy_4s_NumIncome.tif, replace) legend(label(1 "Mean Happiness") label(2 "Linear fit") label(3 "Quadratic fit")) 
graph save FigureA1.tif, replace
	   
drop bins
drop happy_4s_mean

//Figure 1: Happy vs. Log(income)
egen bins=cut(income),at(0(0.2)11)
replace bins=bins+0.1

bysort bins: egen happy_4s_mean=mean(happy_4s)
	  
label variable happy_4s_mean "Mean Happiness"
/* Linear fit and qfit together
twoway (scatter happy_4s_mean bins) (lfit happy_4s income)(qfit happy_4s income), ///
//       ytitle(Mean Happiness) xtitle(Log(income)) saving(MeanHappiness_IncomeBin.tif, replace) legend(label(1 "Mean Happiness") label(2 "Linear fit") label(3 "Quadratic fit")) 
*/
//only linear fit
twoway (scatter happy_4s_mean bins) (lfit happy_4s income), ///
       ytitle(Mean Happiness) xtitle(Log(income)) saving(MeanHappiness_IncomeBin.tif, replace) legend(label(1 "Mean Happiness") label(2 "Linear fit"))  
	   	   
graph save Figure1-lfit.tif, replace
	   
drop bins
drop happy_4s_mean

//Section "Data and Variables": Happiness Distribution
tab happy_4s

//Table 2: Happiness Distribution by Income Brackets
tab happy_4s  income_bracket

//Table 1: Summary Statistics
gen IVincome_num = exp(IVincome)-1

sum happy_4s income income_num IVincome_num age agesquared male BirthOrder_ElderDummy birth_weight VeryEarly_disease_d early_disease_d edu working work_over_40h  married divorced widowed population  
sum happy_4s income income_num IVincome_num age agesquared male BirthOrder_ElderDummy birth_weight VeryEarly_disease_d early_disease_d edu working work_over_40h  married divorced widowed population if mz
sum happy_4s income income_num IVincome_num age agesquared male BirthOrder_ElderDummy birth_weight VeryEarly_disease_d early_disease_d edu working work_over_40h  married divorced widowed population  if dz

//OLS and IV with both twins with complete info; use exactly same sample to estimate OLS
//OLS
//Table S1: specifications 1-6; Table 3 Column 1: specification 6 with full controls
quiet reg happy_4s income age agesquared male $citylist,robust cluster(twinpair)
outreg2 using TableS1_OLS.doc,  nolabel replace coefastr se drop(d_*)

quiet reg happy_4s income age agesquared male BirthOrder_ElderDummy birth_weight VeryEarly_disease_d $citylist,robust cluster(twinpair)  
outreg2 using TableS1_OLS.doc,  nolabel append coefastr se drop(d_*)

quiet reg happy_4s income age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu $citylist,robust cluster(twinpair)
outreg2 using TableS1_OLS.doc,  nolabel append coefastr se drop(d_*)

quiet reg happy_4s income age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu working work_over_40h $citylist,robust cluster(twinpair)  
outreg2 using TableS1_OLS.doc,  nolabel append coefastr se drop(d_*)

quiet reg happy_4s income age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu working work_over_40h married divorced widowed  $citylist,robust cluster(twinpair)
outreg2 using TableS1_OLS.doc,  nolabel append coefastr se drop(d_*)

quiet reg happy_4s income age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu working work_over_40h married divorced widowed population  $citylist,robust cluster(twinpair) //(Table 3, Column 1 OLS)
outreg2 using TableS1_OLS.doc,  nolabel append coefastr se drop(d_*)
outreg2 using Table3_MainResults.doc,  nolabel replace coefastr se drop(d_*)

//Table S2-S3: full results of OLS using numberical income with squared term (Table S2), ordered probit (Table S3)
//Table S2
replace income_num = income_num / 10000
gen income_num_squared = income_num * income_num

quiet reg happy_4s income_num age agesquared male $citylist,robust cluster(twinpair)
outreg2 using TableS2_NumbericIncome.doc,  nolabel replace coefastr se drop(d_*)

quiet reg happy_4s income_num income_num_squared age agesquared male $citylist,robust cluster(twinpair)
outreg2 using TableS2_NumbericIncome.doc,  nolabel append coefastr se drop(d_*)

quiet reg happy_4s income_num income_num_squared age agesquared male BirthOrder_ElderDummy birth_weight VeryEarly_disease_d $citylist,robust cluster(twinpair)
outreg2 using TableS2_NumbericIncome.doc,  nolabel append coefastr se drop(d_*)

quiet reg happy_4s income_num income_num_squared age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu $citylist,robust cluster(twinpair)
outreg2 using TableS2_NumbericIncome.doc,  nolabel append coefastr se drop(d_*)

quiet reg happy_4s income_num income_num_squared age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu working work_over_40h  $citylist,robust cluster(twinpair)
outreg2 using TableS2_NumbericIncome.doc,  nolabel append coefastr se drop(d_*)

quiet reg happy_4s income_num income_num_squared age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu working work_over_40h married divorced widowed $citylist,robust cluster(twinpair)
outreg2 using TableS2_NumbericIncome.doc,  nolabel append coefastr se drop(d_*)
 
quiet reg happy_4s income_num income_num_squared age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu working work_over_40h married divorced widowed population  $citylist,robust cluster(twinpair)
outreg2 using TableS2_NumbericIncome.doc,  nolabel append coefastr se drop(d_*)

//Table S3:
quiet oprobit happy_4s income age agesquared male $citylist,robust cluster(twinpair)
outreg2 using TableS3_OProbit.doc,  nolabel replace coefastr se drop(d_*)

quiet oprobit happy_4s income age agesquared male BirthOrder_ElderDummy birth_weight VeryEarly_disease_d $citylist,robust cluster(twinpair)
outreg2 using TableS3_OProbit.doc,  nolabel append coefastr se drop(d_*)

quiet oprobit happy_4s income age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu $citylist,robust cluster(twinpair)
outreg2 using TableS3_OProbit.doc,  nolabel append coefastr se drop(d_*)

quiet oprobit happy_4s income age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu working work_over_40h  $citylist,robust cluster(twinpair)
outreg2 using TableS3_OProbit.doc,  nolabel append coefastr se drop(d_*)

quiet oprobit happy_4s income age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu working work_over_40h married divorced widowed $citylist,robust cluster(twinpair)
outreg2 using TableS3_OProbit.doc,  nolabel append coefastr se drop(d_*)
 
quiet oprobit happy_4s income age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu working work_over_40h married divorced widowed population  $citylist,robust cluster(twinpair)
outreg2 using TableS3_OProbit.doc,  nolabel append coefastr se drop(d_*)

//Table 3, Column 2: IV
ivreg2 happy_4s (income=IVincome) age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu working work_over_40h married divorced widowed population  $citylist,robust cluster(twinpair) first  
outreg2 using Table3_MainResults.doc,  nolabel append coefastr se drop(d_*)

//no robust and cluster option to allow for hausman test
quiet ivreg2 happy_4s (income=IVincome) age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu working work_over_40h married divorced widowed population  $citylist 
est store IV

//IVFE-1: using xtivreg
xtset twinpair twcode
//MZ
xtivreg happy_4s (income=IVincome) age agesquared male BirthOrder_ElderDummy birth_weight early_disease_d edu working work_over_40h married divorced widowed population  if mz==1,fe  first
est store IVFE_MZ
hausman IVFE_MZ IV,  force

//using first difference
collapse $alllist, by(twinpair) 
		 
***  Table 3, Column 3：IVFE-MZ ***
ivreg2 Dhappy_4s (Dincome=IVDincome)  DBirthOrder_ElderDummy Dbirth_weight Dearly_disease_d Dedu Dworking Dwork_over_40h Dmarried Ddivorced Dwidowed Dpopulation if mz==1,noc robust cluster(twinpair) first 
outreg2 using Table3_MainResults.doc,  nolabel append coefastr se 

***  Table 3, Column 4：IVFE-DZ ****
ivreg2 Dhappy_4s (Dincome=IVDincome)  DBirthOrder_ElderDummy Dbirth_weight Dearly_disease_d Dedu Dworking Dwork_over_40h Dmarried Ddivorced Dwidowed Dpopulation if dz==1,noc robust cluster(twinpair) first 
outreg2 using Table3_MainResults.doc,  nolabel append coefastr se 

******************** Table 4: Robustnesss Checks for IVFE-MZ Results *********
//happiness as 3-scales
ivreg2 Dhappiness (Dincome=IVDincome)  DBirthOrder_ElderDummy Dbirth_weight Dearly_disease_d Dedu Dworking Dwork_over_40h Dmarried Ddivorced Dwidowed Dpopulation if mz==1,noc robust cluster(twinpair)
outreg2 using Table4_Robustnesss.doc,  nolabel replace coefastr se 

//happiness dummy
ivreg2 Dhappy_d (Dincome=IVDincome)  DBirthOrder_ElderDummy Dbirth_weight Dearly_disease_d Dedu Dworking Dwork_over_40h Dmarried Ddivorced Dwidowed Dpopulation if mz==1,noc robust cluster(twinpair) 
outreg2 using Table4_Robustnesss.doc,  nolabel append coefastr se 

//only born before 1979
ivreg2 Dhappy_4s (Dincome=IVDincome)  DBirthOrder_ElderDummy Dbirth_weight Dearly_disease_d Dedu Dworking Dwork_over_40h Dmarried Ddivorced Dwidowed Dpopulation if mz==1 & birth_year<1979,noc robust cluster(twinpair) 
outreg2 using Table4_Robustnesss.doc,  nolabel append coefastr se

//married only
drop if Dedu_spouse==.
//use spouse edu at survey
ivreg2 Dhappy_4s (Dincome =IVDincome) DBirthOrder_ElderDummy Dbirth_weight Dearly_disease_d Dedu Dedu_spouse Dworking Dwork_over_40h Dpopulation if mz, noc robust
outreg2 using Table4_Robustnesss.doc,  nolabel append coefastr se

//use spouse edu at marriage
ivreg2 Dhappy_4s (Dincome =IVDincome) DBirthOrder_ElderDummy Dbirth_weight Dearly_disease_d Dedu Dedu_spouse_atmarriage Dworking Dwork_over_40h Dpopulation if mz, noc robust
//outreg2 using Table4_Robustness_MarriedSample.doc,  nolabel append coefastr sexd	 	

//Table 5: See Heterogeity code file

log close
clear
exit